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Abstract 

The joint probability density function (PDF) of turbulent velocity and concentration of a passive scalar in an urban 
street canyon is computed using a newly developed particle-in-cell Monte Carlo method. Compared to moment 
closures, the PDF methodology provides the full one-point one-time PDF of the underlying fields containing all 
higher moments and correlations. The small-scale mixing of the scalar released from a concentrated source at the 
street level is modelled by the interaction by exchange with the conditional mean (IECM) model, with a micro-mixing 
time scale designed for geometrically complex settings. The boundary layer along no-slip walls (building sides and 
tops) is fully resolved using an elliptic relaxation technique, which captures the high anisotropy and inhomogeneity 
of the Reynolds stress tensor in these regions. A less computationally intensive technique based on wall functions to 
represent boundary layers and its effect on the solution are also explored. The calculated statistics are compared to 
experimental data and large-eddy simulation. The present work can be considered as the first example of computation 
of the full joint PDF of velocity and a transported passive scalar in an urban setting. The methodology proves 
successful in providing high level statistical information on the turbulence and pollutant concentration fields in 
complex urban scenarios. 

Key words: Langevin equation; Monte-Carlo method; Probability density function method; Scalar dispersion; Urban-scale 
turbulence 



1. Introduction 



Regulatory bodies, architects and town planners increasingly use computer models to assess ventilation and 
occurrences of hazardous pollutant concentrations in cities. These models are mostly based on the Reynolds- 
averaged Navier-Stokes (RANS) equations or, more recently, large-eddy simulation (LES) techniques. Both 
of these approaches require a series of modelling assumptions, including most commonly the eddy-viscosity 
and gradient-diffusion hypotheses. The inherent limitations of these approximations, even in the simplest 
engineering flows, are well known, and detailed for example by Pope (2000). Although the effects of the 
assumptions are more pronounced in RANS than in LES models, where they are confined to the smaller 
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(modelled) scales, there is clearly a need to develop higher-order models. In pollutant dispersion modelling 
it is also desirable to predict extreme events such as peak values or probabilities that concentrations will 
exceed a certain threshold. In other words, a fuller statistical description of the concentration field is required 
(Chatwin and Sullivan, 1993; Kristensen, 1994; Wilson, 1995; Pavageau and Schatzmann, 1999). These issues 
have been explored in grid turbulence and in the unobstructed atmosphere, and models capable of predicting 
higher-order statistics have also appeared (Yee et al., 1994; Yee and Wilson, 2000; Luhar et al., 2000; Cassiani 
and Giostra, 2002; Franzese, 2003; Cassiani et al., 2005a, b), but more research is necessary to extend these 
capabilities to cases of built-up areas. 

Probability density function (PDF) methods have been developed mainly within the combustion engineer- 
ing community as an alternative to moment closure techniques to simulate chemically reactive turbulent flows 
(Lundgren, 1969; Pope, 1985; Dopazo, 1994). Because many-species chemistry is high-dimensional and highly 
non-linear, the biggest challenge in reactive flows is to adequately model the chemical source term. In PDF 
methods, the closure problem is raised to a statistically higher level by solving for the full PDF of the turbu- 
lent flow variables instead of its moments. This has several benefits: convection, the effect of mean pressure, 
viscous diffusion and chemical reactions appear in closed form in the PDF transport equation. Therefore 
these processes are treated mathematically exactly without closure assumptions, eliminating the need for 
gradient-transfer approximations. The effects of fluctuating pressure, dissipation of turbulent kinetic energy 
and small-scale mixing of scalars still have to be modelled. The rationale is that, since the most important 
physical processes are treated exactly, the errors introduced by modelling assumptions for less important 
processes amount to a smaller departure from reality. Moreover, the higher level description provides more 
information that can be used in the construction of closure models. 

The PDF transport equation is a high-dimensional scalar equation. Although techniques of solution 
based on stochastic Eulerian methods have been developed (Valiho, 1998; Mobus et al., 2001; Soulard 
and Sabel'nikov, 2006), Lagrangian Monte Carlo methods arc a more natural choice because their compu- 
tational cost increases only linearly with the dimensionality of the problem, favourably comparing to the 
more traditional finite difference, finite volume or finite element methods. 

The numerical development of Lagrangian PDF methods has mainly centred around three distinctive 
approaches, all of them representing a finite ensemble of fluid particles with Lagrangian particles. A common 
approach is the stand-alone Lagrangian method, where the flow is represented by particles whereas the 
Eulerian statistics are obtained using kernel estimation (Pope, 2000; Fox, 2003). Another technique is the 
hybrid methodology, which builds on existing Eulerian computational fluid dynamics (CFD) codes based 
on moment closures (Muradoglu et al., 1999, 2001; Jenny et al., 2001; Rembold and Jenny, 2006). Hybrid 
methods use particles to solve for certain quantities and provide closures for the Eulerian moment equations 
using the particle/PDF methodology. A more recent approach is the self-consistent non-hybrid method 
(Bakosi et al., 2007, 2008), which also employs particles to represent the flow, and uses the Eulerian grid only 
to solve for inherently Eulerian quantities (such as the mean pressure) and for efficient particle tracking. Since 
the latter two approaches extensively employ Eulerian grids, they are particle-in-ccll methods (Grigoryev 
et al., 2002). 

The current study presents an application of the non-hybrid method to a simplified urban-scale case where 
pollution released from a concentrated line source between idealized buildings is simulated and results are 
compared to data from wind-tunnel experiments and LES. 

PDF methods in atmospheric modelling mostly focus on the simulation of passive pollutants, wherein 
the velocity field (mean and turbulence) is either assumed or obtained from experiments (Sawford, 2004, 
2006; Cassiani et al., 2005a, b, 2007). In contrast, the current model directly computes the joint PDF of the 
turbulent velocity, characteristic turbulent frequency and scalar concentration, extending the use of PDF 
methods in atmospheric modelling to represent more physics at a higher statistical level. A computed full 
joint PDF also has the advantage of providing information on the uncertainty originating from turbulence 
on a physically sound basis. 

In this study the turbulent boundary layers developing along solid walls are treated in two different ways: 
either fully resolved or via the application of wall functions (i.e. the logarithmic "law of the wall"). The full 
resolution is obtained using Durbin's elliptic relaxation technique (Durbin, 1993), which was incorporated 
into the PDF methodology by Dreeben and Pope (1997a, 1998). This technique allows for an adequate 
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representation of the near- wall low- Reynolds-number effects, such as the high inhomogeneity and anisotropy 
of the Reynolds stress tensor and wall-blocking. Wall conditions for particles based on the logarithmic "law 
of the wall" in the PDF framework have also been developed by Dreeben and Pope (1997b). These two 
types of wall treatments are examined in terms of trade-off between computational cost and performance, 
addressing the question of how important it is to adequately resolve the boundary layers along solid walls 
in order to obtain reasonable scalar statistics. 

At the urban scale the simplest settings to study turbulent flow and dispersion patterns are street canyons. 
Due to increasing concerns for environmental issues and air quality standards in cities, a wide variety of 
canyon configurations and release scenarios have been studied both experimentally (Hoydysh et al., 1974; 
Wedding et al., 1977; Rafailids and Schatzmann, 1995; Meroney et al., 1996; Pavageau and Schatzmann, 
1999) and numerically (Lee and Park, 1994; Johnson and Hunter, 1995; Baik and Kim, 1999; Huang et al., 
2000; Liu and Barth, 2002). Street canyons have a simple flow geometry, they can be studied in two dimen- 
sions, and a wealth of experimental and modelling data are available for different street-width to building- 
height ratios. This makes them ideal candidates for testing a new urban pollution dispersion model. We 
validate the computed velocity and scalar statistics with the LES results of Liu and Barth (2002) and the 
wind-tunnel measurements of Meroney et al. (1996); Pavageau (1996); Pavageau and Schatzmann (1999). 
The experiments have been performed in the atmospheric wind tunnel of the University of Hamburg, where 
the statistics of the pollutant concentration field have been measured in an unusually high number of loca- 
tions in order to provide fine details inside the street canyon. 

The rest of the paper is organized as follows. In Section 2, the exact and modelled governing equations 
are presented. Several statistics are compared to experimental data and large-eddy simulations in Section 3. 
Finally, Section 4 draws some conclusions and elaborates on possible future directions. 



2. Governing equations 



We write the Eulerian governing equations for a passive scalar released in a viscous, Newtonian, incom 
pressible flow as 

dx t 



0, 



H¥i + u + 1— - v v 2 u 

dt 3 dxj p dxi 



dt 



Ui 



dxj 



TV 2 



(1) 
(2) 
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where Ui, P, p, v, <p an d T are the Eulerian velocity, pressure, constant density, kinematic viscosity, scalar 
concentration and scalar diffusivity, respectively. Based on this system of equations the exact transport 
equation that governs the one-point, one-time Eulerian joint PDF of velocity and concentration f(V, x, i) 
can be written as (Pope, 1985, 2000), 
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d 2 f 1 d{P) df 
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(4) 



where V and ip denote the sample space variables of the stochastic velocity U(x,t) and concentration 
4>{x,i) fields, respectively, and the pressure P is decomposed into its mean (P) and fluctuation part p. In 
Equation (4) the physical processes of advection (second term on the left), viscous diffusion (first term on 
the right) and transport of / in velocity space by the mean pressure gradient (second term on the right) 
are represented mathematically exactly. The last three terms in the form of conditional expectations have 
to be modelled: these are respectively the effect of dissipation of turbulent kinetic energy, the effect of 
fluctuating pressure and the small-scale diffusion of the scalar. After appropriate modelling of the unclosed 
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terms, Equation (4) can, in principle, be solved with a traditional numerical method. However, the high- 
dimensionality makes Monte Carlo methods more appealing. In particular, because Equation (4) is a Fokkcr- 
Planck equation, it can be written as an equivalent system of stochastic differential equations (SDEs) (van 
Kampen, 2004). We use the generalized Langcvin model (GLM) of Haworth and Pope (1986) for the velocity 
increment, and the interaction by exchange with the conditional mean (IECM) model for the scalar. The 
physics and characteristics of the IECM model are discussed in detail by Fox (1996); Pope (1998) and 
Sawford (2004). Thus our system of SDEs that solve Equation (4) is written as 

dXi =Uidt + (2vf 12 dW h (5) 

dUi = ~li£^ dt + 2i/ S£ d< + {2v)1 ' 2 ^ dWj + Gij {Uj - {Uj)) At + {cq£)1/2 dw *> (6) 

dip = --L^-(cf>\U = V))dt. (7) 

Equation (5) governs the Lagrangian particle position A^, and it consists of advection by the instantaneous 
particle velocity lAi and molecular diffusion represented by the isotropic Wiener process dWi, which is a 
given Gaussian process with zero mean and variance dt. The particle velocity Ui is governed by Equation (6). 
The second and third terms on the right-hand side are a direct consequence of the particular Lagrangian 
representation of the viscous diffusion in Equation (5) (for a derivation, see Dreeben and Pope, 1997a), thus 
the first three terms on the right-hand side govern the particle velocity in a laminar flow mathematically 
exactly. The last two terms involving the tensor Gy and Co jointly model the effect of pressure redistribution 
and anisotropic dissipation of turbulent kinetic energy. Gij is a second-order tensor function of the mean 
velocity gradients d(Ui) /dxj, the Reynolds stresses (uiUj) and the rate of dissipation of turbulent kinetic 
energy e, while Co is a positive constant. Note that the Wiener process dWi appearing in both Equations (5) 
and (6) is the same process, i.e. the same exact series of random numbers, and is independent of the other 
process dW(. The concentration of the passive scalar is governed by Equation (7), which represents the 
physical process of diffusion by relaxation of the particle scalar ip towards the velocity-conditioned scalar 
mean (<fi\U = V) = (<j>\V) with a timescale i m . Note that the Eulerian statistics denoted by angled brackets 
(•} are to be evaluated at the fixed particle locations Xi, and in particle-in-cell methods, this is usually 
achieved using an Eulerian grid and computing ensemble averages in grid elements and/or around vertices. 
The energy dissipation rate is defined as 

e = {lo) (k + vC^{uj)) , (8) 

where Cr is a model constant. We adopt the model of van Slooten ct al. (1998) for the stochastic characteristic 
turbulent frequency oj 

doj = -C a (w)(w - {u))dt - S u (u)ujdt + (2C 3 C 4 (Lj) 2 Lj) 1/2 dW, (9) 
where is a source/sink term for the mean turbulent frequency 
V 

Su> = C W 2 — C w i — , (10) 

where V = —(uiUj)d(Ui) / 'dxj is the production of turbulent kinetic energy, dW is a scalar valued Wiener 
process, while CjjC^Cji and C U 2 arc model constants. To define the micro-mixing time scale for a scalar 
released from a concentrated source in a geometrically complex flow domain we follow Bakosi et al. (2007, 
2008) and specify the inhomogeneous t m as 



t m (x) 



(ii) 



where ro denotes the radius of the source, U c is a characteristic velocity at x that we take as the Euclidean 
norm of the mean velocity vector, Xq is the location of the source, and C s and Ct are model constants. 
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2.1. Elliptic relaxation modelling of near-wall turbulence 



The turbulence model GLM (represented by the last two terms of Equation (6) for the velocity increments) 
is in fact a family of models. A specific definition of Gij and Cq corresponds to a particular closure from 
a wealth of models, many of which can be made equivalent (at the level of second moments) to popular 
Reynolds-stress closures (Pope, 1994). To be able to capture the near- wall low- Reynolds- number effects on 
the Reynolds stresses in fully resolved boundary layers, we follow Durbin (1993) and Dreeben and Pope 
(1998) and specify Gij and Cq through the tensor pij 



S(f= V^ . (12) 



2pjj - e5 lJ 

Hi ~ 2k 

2p ij (u i u j ) 

° Q - 3fc^' (13) 
where k = ^(uiUi) denotes the turbulent kinetic energy and pij is obtained by solving the elliptic equation 

Pij - L 2 W 2 pij = |(1 - C^k^Sij + kHijJ^-, (14) 
where the fourth-order tensor Hijki is given by 

Hijki = (C 2 A V + \^)8 ik 5ji - \^8 u 5 ]k + ~fsb%k&ji ~ Ishi&jk, (15) 

with 



A v 



1; a(ffc) 8 det («,•«,•) , (16) 



the Reynolds stress anisotropy 

and Ci, C2, 75, C„ are model constants. The characteristic length scale L in Equation (14) is defined by the 
maximum of the turbulent and Kolmogorov length scales 

L = C L max \c ( (fc 3/2 /e) ; C v {v 3 /e) 1/4 j , (18) 

with = 1 + l.Snini, where Cl and are model constants, and m is the unit wall- normal of the closest 
wall element pointing outward of the flow domain. For the applied wall-boundary conditions in this fully- 
resolved case the reader is referred to Dreeben and Pope (1998) and Bakosi et al. (2008). More details on the 
inflow and outflow conditions for the mean pressure and on the wall conditions for the tensor pij are given in 
(Bakosi et al., 2008). Equation (14) was developed in conjunction with turbulent channel flow (Durbin, 1993; 
Dreeben and Pope, 1998). Modifications and different forms of this idea have also been proposed (Whizman 
et al., 1996; Dreeben and Pope, 1997a, 1998; Waclawczyk et al., 2004), and an application in channel flow 
with the current non- hybrid model is presented by Bakosi et al. (2007). The model to compute the joint 
PDF of velocity, scalar, and characteristic turbulent frequency is now complete. 



2.2. Wall functions modelling of near-wall turbulence 



When the wall region has sufficient resolution, Gij and Co as defined by Equations (12) and (13) enable the 
model to adequately capture the near-wall effects on the higher-order turbulence statistics. Full resolution 
at the wall is strictly required in certain cases such as, for example, computations of heat transfer at walls 
embedded in a flow or detaching boundary layers with high adverse pressure gradients. In many other realistic 
simulations, however, full resolution of high-Reynolds-number boundary layers is not always possible and 
may not be necessary when the flow details close to walls are not important because the analysis focuses on 
the boundary-layer effects at farther distances. In these cases an alternative option is to model the near-wall 
turbulence using wall functions instead of the elliptic relaxation technique. Employing wall functions for 
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no-slip walls provides a trade-off between the accuracy of fully resolved boundary layers and computational 
speed. Wall functions are widely applied in atmospheric simulations, where full wall resolution is usually 
prohibitively expensive even at the microscale or urban scale (Bacon et al., 2000; Lien et al., 2004). It is 
worth emphasizing that one of the main assumptions used in the development of wall functions is that the 
boundary layer remains attached, which is not always the case in simulations of complex flows. However, 
since wall functions are the only choice for realistic atmospheric simulations, they are still routinely employed 
with reasonable success. 

To investigate the gain in performance and the effect on the results, we implemented the wall treatment 
for complex flow geometries that has been developed for the PDF method by Dreeben and Pope (1997b). 
In this case, the tensor Gy is defined by the simplified Langevin model (SLM) (Haworth and Pope, 1986) 
and Co is simply a constant: 

Gij = - (§ + f Co) (cj>% (19) 

with Co = 3.5. In line with the purpose of wall functions, boundary conditions have to be imposed on 
particles that hit the wall so that their combined effect on the statistics at the first grid point from the 
wall will be consistent with the universal logarithmic wall function in equilibrium flows, i.e. in boundary 
layers with no significant adverse pressure gradients. The development of boundary conditions based on 
wall functions rely on the self-similarity of attached boundary layers close to walls. These conditions are 
applied usually at the first grid point from the wall based on the assumption of constant or linear stress 
distribution. This results in the well-known self-similar logarithmic profile for the mean velocity. For the 
sake of completeness the conditions on the particles developed by Dreeben and Pope (1997b) are repeated 
here. The condition for the wall-normal component of the particle velocity reads 

V R = -V t , (20) 

where the subscripts R and / denote reflected and incident particle properties, respectively. The reflected 
strcamwise particle velocity is given by 

U R = U I + aV I , (21) 

where the coefficient a is determined by imposing consistency with the logarithmic law at the distance of 
the first grid point from the wall, y p : 

= W ' ( } 

where u p is a characteristic velocity scale of the turbulence intensity in the vicinity of y p , defined as 

u P = C]!%y\ (23) 

with C p = 0.09. (U) , (v 2 ) p and k p are, respectively, the mean streamwise velocity, the wall-normal com- 
ponent of the Reynolds stress tensor and the turbulent kinetic energy at y p . In Equation (22) U e is the 
magnitude of the equilibrium value of the mean velocity at y p and is specified by the logarithmic law 

L/ e = ^log(^), (24) 

where k = 0.41 is the Karman constant and the surface roughness parameter E = 8.5 for a smooth wall. 
The friction velocity it* is computed from local statistics as 



u* = \\u% + 7- 



y P d(P) 



7 T = max 



p dx 
0; sign ( (uv) 



with 

d(P) 



dx 



(25) 



(26) 



In Equations (20-25) the streamwise x and wall-normal y coordinate directions are defined according to 
the local tangential and normal coordinate directions of the particular wall element in question. In other 
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Fig. 1. Geometry and Eulerian mesh (consisting of approximately 12,000 cells) for the computation of turbulent street canyon 
with full resolution of the wall boundary layers using elliptic relaxation. The grid is generated by the general purpose mesh 
generator Gmsh (Geuzaine and Remacle, 2009). The positions labelled by bold numbers indicate the sampling locations for 
the passive scalar, equivalent with the combined set of measurement tapping holes of Meroney et al. (1996), Pavageau (1996) 
and Pavageau and Schatzmann (1999). In the zoomed area the refinement is depicted, which ensures an adequate resolution of 
the boundary layer and the vortices forming in the corner. 

words, if the wall is not aligned with the flow coordinate system then the vectors Ui and d(P)/dxi, and the 
Reynolds stress tensor (uiUj), need to be appropriately transformed into the wall-element coordinate system 
before being employed in the above equations. 

The condition on the turbulent frequency is given by 



which completes the description of the wall function approach. 

In summary, the flow is represented by a large number of Lagrangian particles whose position A^, velocity 
Ui, scalar concentration ip and characteristic turbulent frequency ui are governed by Equations (5), (6), (7) 
and (9), respectively. These equations are discretised and advanced in time by the explicit forward Euler- 
Maruyama method (Kloeden and Platen, 1999). The mean pressure, required in Equation (6), is obtained via 
a pressure projection scheme (Bakosi et al., 2008). Full wall resolution is obtained through Equations (12-18), 
while wall functions are applied through Equations (19-28). The pressure-Poisson and elliptic relaxation (14) 
equations are solved using an unstructured Eulerian grid with the finite element method. The grid is also used 
to track particles throughout the domain and to estimate Eulerian statistics using ensemble averaging. In 
practical simulations using PDF methods a few hundred particles per clement is usually employed. Adequate 
stability can already be achieved using as little as 50-100 particles, however, 300-500 particles per elements 
are recommended to exploit the bin structure to compute {4>\V) (Bakosi et al., 2008) and to decrease the 
statistical error. The numerical algorithm and performance issues are detailed in Bakosi ct al. (2008). 

3. Modelling the street canyon 




(27) 



with 



2 



(28) 



5 + f Co + C3 + C U 2 — Cui 



Street canyons are often used to study flow and pollutant dispersion patterns in urban areas. A wealth 
of experimental data for this simplified urban-scale setting is available from wind-tunnel and LES data 
(Meroney et al., 1996; Pavageau and Schatzmann, 1999; Liu and Barth, 2002), making it a natural choice 
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Fig. 2. Geometry and Eulerian mesh (consisting of approximately 500 cells) for the computation of flow in a turbulent street 
canyon with wall functions at Re m 12000. The domain is stripped at no-slip walls so that it does not include the close vicinity 
of the wall at y+ < 30. The positions for sampling the scalar concentrations are the same as in Figure 1. 

to validate the current, newly developed method. We will simulate the "urban roughness" case of Meroney 
et al. (1996), which is a model for a series of street canyons in the streamwise direction. The simulations 
are performed for statistically two-dimensional flow geometry, with periodic inflow and outflow boundary 
conditions in the free stream above the buildings (i.e. the particles crossing the outflow boundary are re- 
injected at the inflow boundary). The Reynolds number based on the maximum free stream velocity Uq and 
the building height H was Re rss 12000. This corresponds to Re T rj 600 based on the friction velocity and 
the free stream height, h = H/2, if the free stream above the buildings is considered as the lower part of an 
approximate fully developed turbulent channel flow. The velocity-conditioned scalar mean (</>|V) required 
in Equation (7) has been computed using the general method described by Bakosi ct al. (2007) using a bin 
structure of (5 x 5 x 5). After the flow has reached a statistically stationary state, time averaging is used 
to collect velocity statistics and a continuous scalar is released from a street level line source at the centre 
of the canyon (corresponding to a point source in two dimensions). Particles travelling through the source 
of diameter 0.01//i are assigned a unit source strength. The scalar field is also time averaged after it has 
reached a stationary state. 

The simulations with the full resolution model have been made with the constants given in Table 2, using 

Table 1 



Concentration sampling locations at building walls and tops according to the experimental measurement holes of Meroney et al. 
(1996), Pavagcau and Schatzmann (1999) and Pavageau (1996). See also Figure 1. 
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Table 2 

Constants for modelling 


the joint PDF of velocity, characteristic 


turbulent frequency and transported passive scalar. 
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Fig. 3. Velocity vectors (first row) and iso-contours of turbulent kinetic energy (second row) of the fully developed turbulent 
street canyon at Re £S 12000 based on the maximum free stream velocity Uq and the building height H. Left - full resolution 
with elliptic relaxation, right - coarse simulation with wall functions. 

300 particles per element. The Eulerian mesh used for this simulation is displayed in Figure 1, which shows 
the considerable refinement along the building walls and tops necessary to resolve the boundary layers. In 
this case, the high anisotropy and inhomogeneity of the Reynolds stress tensor in the vicinity of walls are 
captured by the elliptic relaxation technique, using Equation (14). 

The simulations using wall functions were performed on the Eulerian mesh displayed in Figure 2, also 
using 300 particles per element. The particle-boundary conditions described in Sec. 2.2 were implemented 
for arbitrary geometry. Note that the first grid point where the boundary conditions based on wall functions 
are to be applied should not be closer to the wall than y + = yu T jv = 30, where y + is the non-dimensional 
distance from the wall in wall units, but should still be sufficiently close to the wall to lie in the inertial 
sublayer (Drecben and Pope, 1997b). Accordingly, the grid in Figure 2 only contains the domain stripped 
from the wall region at y + < 30. 

Turbulence and scalar statistics are obtained entirely from the particles that represent both the flow itself 
and the scalar concentration field. The Eulerian meshes displayed in Figure 1 for the full resolution and in 
Figure 2 for the wall function cases are used to extract the statistics, to track the particles throughout the 
domain and to solve the Eulerian equations, namely Equation (14) and the mean-pressurc-Poisson equation 
in the fully resolved case and only the latter in the wall function case. 

In Figure 3, the mean velocity vector field and the iso-contours of the turbulent kinetic energy are displayed 
for both fully resolved and wall function simulations. It is apparent that the full resolution captures even the 
smaller counter-rotating eddies at the internal corners of the canyon, while the coarse grid resolution with 
wall functions only captures the overall flow pattern characteristic of the flow, such as the large steadily 
rotating eddy inside the canyon. The turbulent kinetic energy field is captured in a similar manner. Both 
methods reproduce the highest turbulence activity at the building height above the canyon, with a maximum 
at the windward building corner. The full resolution simulation shows a more detailed spatial distribution 
of energy, whereas the coarse resolution of the wall-function simulation still allows tha capture of the overall 
pattern. 
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Fig. 4. Dimensionless turbulent intensities (u 2 ) ' /Co (first column) and ("if 2 ^ ' /Uq (second column) computed using full 
wall resolution (first row) and using wall functions (second row) at Re ~ 12000 compared with the LES results (third row) of 
Liu and Barth (2002). 

In Figure 4 the normalized turbulent intensities (u 2 ) 1/2 /Uq and (w 2 ) 1/2 /Uq are displayed for both simula- 
tion cases and compared with the LES results of Liu and Barth (2002). In the large eddy simulations the 
filtered momentum equations are solved by the Galerkin finite clement method using brick three-dimensional 
elements, while the residual stresses are modelled using the Smagorinsky closure. 

The full resolution simulation shows very good agreement with LES. The contour plots of (u 2 ) 1 ' 2 /Uq 
correctly display two local maxima, at the windward external and at the leeward internal corners. The 
contour plots of (w 2 ) 1/2 /Uq show distributed high values at the building level above the canyon, along the 
windward internal corner and wall, and at the street level downstream of the source. By contrast, the wall- 
function contour plots are in general less detailed, failing to reproduce the internal maximum of [u 2 ) 1/2 /Uq, 
and showing a more uniform representation of {w 2 ) 1/2 /Uq. 

Several wind-tunnel measurements have been carried out for a scalar released from a street level continuous 
line source at the centre of the canyon, providing concentration statistics above the buildings, at the walls, 
and inside the canyon (Meroncy ct al., 1996; Pavageau, 1996; Pavageau and Schatzmann, 1999). To examine 
the concentration values along the building walls and tops, we sampled the computed mean concentration 
field at the locations depicted in Figure 1 and listed in Table 1. 

The excellent agreement of the results using both full resolution and wall functions with a number of 
experiments is shown in Figure 5. The concentration peak is precisely captured at the internal leeward 
corner and the model accurately reproduces the pattern of concentration along both walls including the 
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Fig. 5. Distribution of mean concentrations at the boundary of the street canyon. The experimental data arc in terms of the 
ratio CU re {H L/Q B , where C is the actual measured mean concentration (ppm). U IC { is the free-stream mean velocity (ms _1 ) 
taken at the reference height y rc f ?» 11H and Q a /L is the line source strength (m 2 s _1 ) in which Q s denotes the scalar flow 
rate and L is the source length. The calculation results are scaled to the concentration range of the experiments. References 
for experimental data: a Meroney ct al. (1996); o, v, Pavageau and Schatzmann (1999); □ Pavageau (1996). See also Figure 1 
and Table 1 for the measurement locations. 

higher values along the leeward wall. 

In Figure 6, the first two statistical moments of the concentration inside the canyon are compared with 
experimental data and LES. The agreement with observations indicates that both the fluid dynamics and 
the micro-mixing components of the model provide a good representation of the real field. 

Because the one-point one-time joint PDF contains all higher-order statistics and correlations of the 
velocity and scalar fields resulting from a close, low-level interaction between the two fields, a great wealth 
of statistical information is available for atmospheric transport and dispersion calculations. As an example, 
in Figure 7 the cumulative distribution functions (CDF) of scalar concentration fluctuations 



are depicted after time averaging at selected locations of the domain for the full resolution case. No experi- 
mental data are available to assess the distributions in Figure 7, although the irregular shape of the CDF at 
x = 3, y = 2, which corresponds to a multimodal PDF, is likely to be an artifact of the micro-mixing model. 

The performance gain obtained by applying wall functions as opposed to full resolution was about two 
orders of magnitude already at our moderate Reynolds number. The gain for higher Reynolds numbers is 
expected to increase faster than linearly. 

4. Discussion 

We have used an Eulcrian unstructured grid, consisting of triangular clement types, to estimate Eulcrian 
statistics, to track particles throughout the domain, and to solve for inherently Eulcrian quantities in con- 
junction with a PDF method. The boundary layers developing close to solid walls are fully captured with an 
elliptic relaxation technique, but can also be represented by wall functions, which use a coarser grid resolu- 
tion and require significantly less particles, resulting in substantial savings in computational cost. We found 
that the one-point statistics of the joint PDF of velocity and scalar are well-captured by the wall function 




(29) 
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Fig. 6. Comparison of the spatial distribution of the normalized mean CU IC {H L/Q B (left column) and variance 
c 2 ^(U Ie fHL/Q B ) 2 (right column) of the scalar released at the centre of the street level. The normalization and the scal- 
ing of the calculated results are the same as in Figure 5. First row - PDF calculations with full wall resolution, second row 
- PDF calculations with wall functions, third row — experimental data of Pavageau and Schatzmann (1999) and fourth row - 
LES calculations of Liu and Barth (2002). 



approximation. In view of its affordable computational load and reasonable accuracy, this approximation 
appears to hold a realistic potential for application of the PDF method in atmospheric simulations, where 
the natural extension of the work is the implementation of the model in three spatial dimensions. 

In hybrid PDF models developed for complex chemically reacting flows, numerical treatments for boundary 
conditions have been included for symmetric, inflow, outflow and free-slip walls employing the ghost-cell 
approach common in finite volume methods (Rembold and Jenny, 2006). The representation of no-slip 
boundaries adds a significant challenge to the above cases. This is partly due to the increased computational 
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Fig. 7. Cumulative distribution functions of scalar concentration fluctuations (left) at x = 3, y = 0.2 and (right) at x = 3, 
y = 2 using full resolution at walls. 

expense entailed by the higher Eulerian grid resolution that is required to fully resolve the boundary layers. 
In addition, there is an increased complexity in specifying the no-slip particle conditions for both fully 
resolved and wall function representations. We presented an implementation of both approaches to treat no- 
slip boundaries with unstructured grids in conjunction with the finite clement method. This also obviates 
further complications with ghost cells. 

In the case of full wall resolution we employed the Lagrangian equivalent of a modified isotropisation of 
production (IP) model as originally suggested by Dreeben and Pope (1998). The elliptic relaxation technique, 
however, allows for the application of any turbulence model developed for high-Reynolds-number turbulence 
(Durbin, 1993; Whizman et al., 1996). The standard test case for developing near- wall models is the fully 
developed turbulent channel flow. In this case, we explored the simpler Rotta (1951) model, which is the 
Eulerian equivalent of the simplified Langevin model (SLM) in the Lagrangian framework (Pope, 1994). 
This is simply achieved by eliminating the term involving the fourth-order tensor Hijki from the right-hand 
side of Equation (14). While the SLM makes no attempt to represent the effect of rapid pressure (Pope, 
2000) (in fact it is strictly correct only in decaying homogeneous turbulence), it is widely applied due to its 
is simplicity and robustness. Our experience showed a slight degradation of the computed velocity statistics 
(as compared to direct numerical simulation) using SLM for the case of channel flow. Since we experienced 
no significant increase in computational expense or decrease in numerical stability, we retained the original 
IP model. 

Similarly, in the case of wall functions, several choices are available regarding the employed turbulence 
model. The methodology developed by Dreeben and Pope (1997b) uses the SLM, but it is general enough 
to include other more complex closures, such as the Haworth and Pope (1986, 1987) models (HP1 and 
HP2), the different variants of the IP models (IPMa, IPMb, LIPM) (Pope, 1994) or the Lagrangian version 
of the SSG model of Speziale et al. (1991). All these closures can be collected under the umbrella of the 
generalized Langevin model, by specifying its constants as described by Pope (1994). These models have been 
all developed for high-Rcynolds-number turbulence and need to be modified in the vicinity of no-slip walls. 
Including them in the wall-function formulation is possible by specifying the reflected particle frequency at 
the wall as 



instead of Equation (28). This involves the additional computation of the statistics (uv) and {u>v 2 } at y P7 
which does not increase the computational cost significantly, but may result in a numerically less stable 
condition since the (originally constant) parameter (3 that appears using the SLM has been changed to 
a variable that fluctuates during simulation. We implemented and tested all the above turbulence models 
using the wall function technique. Without any modification of the model constants we found the IPMa and 
SLM to be the most stable, providing very similar results. Thus we retained the original (and simplest) SLM 



U) R = cxp(-2Vi(ujv) p /(ujv 2 ) p ) 



(30) 
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along with Equation (28). 

The most widely employed closure to model the small scale mixing of the passive scalar in the Lagrangian 
framework is the interaction by exchange with the mean (IEM) model of Villcrmaux and Dcvillon (1972) and 
Dopazo and O'Brien (1974). This simple and efficient model, however, fails to comply with several physical 
constraints and desirable properties of an ideal mixing model (Fox, 2003). The interaction by exchange 
with the conditional mean (IECM) model overcomes some of the difficulties inherent in the IEM model. 
In this study we justify the use of the IECM model by its being more physical and more accurate, but we 
acknowledge that it markedly increases the computational cost (Bakosi et al., 2008). 
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